ANCOVA in Vocabulary taught (Vocabulary taught)

Geiser C. Challco geiser@alumni.usp.br

NOTE:

Setting Initial Variables

dv = "score.vocab.ensinado"
dv.pos = "score.vocab.ensinado.pos"
dv.pre = "score.vocab.ensinado.pre"

fatores2 <- c("genero","zona.participante","zona.escola","score.vocab.ensinado.quintile")
lfatores2 <- as.list(fatores2)
names(lfatores2) <- fatores2

fatores1 <- c("grupo", fatores2)
lfatores1 <- as.list(fatores1)
names(lfatores1) <- fatores1

lfatores <- c(lfatores1)

color <- list()
color[["prepost"]] = c("#ffee65","#f28e2B")
color[["grupo"]] = c("#bcbd22","#fd7f6f")
color[["genero"]] = c("#FF007F","#4D4DFF")
color[["zona.escola"]] = c("#AA00FF","#00CCCC")
color[["zona.participante"]] = c("#AA00FF","#00CCCC")

level <- list()
level[["grupo"]] = c("Controle","Experimental")
level[["genero"]] = c("F","M")
level[["zona.escola"]] = c("Rural","Urbana")
level[["zona.participante"]] = c("Rural","Urbana")

# ..

ymin <- 0
ymax <- 0

ymin.ci <- 0
ymax.ci <- 0


color[["grupo:genero"]] = c(
  "Controle:F"="#ff99cb", "Controle:M"="#b7b7ff",
  "Experimental:F"="#FF007F", "Experimental:M"="#4D4DFF",
  "Controle.F"="#ff99cb", "Controle.M"="#b7b7ff",
  "Experimental.F"="#FF007F", "Experimental.M"="#4D4DFF"
)
color[["grupo:zona.escola"]] = c(
  "Controle:Rural"="#b2efef","Controle:Urbana"="#e5b2ff",
  "Experimental:Rural"="#00CCCC", "Experimental:Urbana"="#AA00FF",
  "Controle.Rural"="#b2efef","Controle.Urbana"="#e5b2ff",
  "Experimental.Rural"="#00CCCC", "Experimental.Urbana"="#AA00FF"
)
color[["grupo:zona.participante"]] = c(
  "Controle:Rural"="#b2efef","Controle:Urbana"="#e5b2ff",
  "Experimental:Rural"="#00CCCC", "Experimental:Urbana"="#AA00FF",
  "Controle.Rural"="#b2efef","Controle.Urbana"="#e5b2ff",
  "Experimental.Rural"="#00CCCC", "Experimental.Urbana"="#AA00FF"
)

for (coln in c(
  "palavras.lidas","score.compreensao","tri.compreensao",
  "score.vocab","tri.vocab",
  "score.vocab.ensinado","tri.vocab.ensinado","score.vocab.nao.ensinado","tri.vocab.nao.ensinado",
  "score.CLPP","tri.CLPP","score.CR","tri.CR",
  "score.CI","tri.CI","score.TV","tri.TV","score.TF","tri.TF","score.TO","tri.TO")) {
  color[[paste0(coln,".quintile")]] = c("#BF0040","#FF0000","#800080","#0000FF","#4000BF")
  level[[paste0(coln,".quintile")]] = c("1st quintile","2nd quintile","3rd quintile","4th quintile","5th quintile")
  color[[paste0("grupo:",coln,".quintile")]] = c(
    "Experimental.1st quintile"="#BF0040", "Controle.1st quintile"="#d8668c",
    "Experimental.2nd quintile"="#FF0000", "Controle.2nd quintile"="#ff7f7f",
    "Experimental.3rd quintile"="#8fce00", "Controle.3rd quintile"="#ddf0b2",
    "Experimental.4th quintile"="#0000FF", "Controle.4th quintile"="#b2b2ff",
    "Experimental.5th quintile"="#4000BF", "Controle.5th quintile"="#b299e5",
    
    "Experimental:1st quintile"="#BF0040", "Controle:1st quintile"="#d8668c",
    "Experimental:2nd quintile"="#FF0000", "Controle:2nd quintile"="#ff7f7f",
    "Experimental:3rd quintile"="#8fce00", "Controle:3rd quintile"="#ddf0b2",
    "Experimental:4th quintile"="#0000FF", "Controle:4th quintile"="#b2b2ff",
    "Experimental:5th quintile"="#4000BF", "Controle:5th quintile"="#b299e5")
}


gdat <- read_excel("../data/data.xlsx", sheet = "vocabulario.st")



dat <- gdat
dat$grupo <- factor(dat[["grupo"]], level[["grupo"]])
for (coln in c(names(lfatores))) {
  dat[[coln]] <- factor(dat[[coln]], level[[coln]][level[[coln]] %in% unique(dat[[coln]])])
}
dat <- dat[which(!is.na(dat[[dv.pre]]) & !is.na(dat[[dv.pos]])),]
dat <- dat[,c("id",names(lfatores),dv.pre,dv.pos)]

dat.long <- rbind(dat, dat)
dat.long$time <- c(rep("pre", nrow(dat)), rep("pos", nrow(dat)))
dat.long$time <- factor(dat.long$time, c("pre","pos"))
dat.long[[dv]] <- c(dat[[dv.pre]], dat[[dv.pos]])


for (f in c("grupo", names(lfatores))) {
  if (is.null(color[[f]]) && length(unique(dat[[f]])) > 0) 
      color[[f]] <- distinctColorPalette(length(unique(dat[[f]])))
}
for (f in c(fatores2)) {
  if (is.null(color[[paste0("grupo:",f)]]) && length(unique(dat[[f]])) > 0)
    color[[paste0("grupo:",f)]] <- distinctColorPalette(length(unique(dat[["grupo"]]))*length(unique(dat[[f]])))
}

ldat <- list()
laov <- list()
lpwc <- list()
lemms <- list()

Descriptive Statistics of Initial Data

df <- get.descriptives(dat, c(dv.pre, dv.pos), c("grupo"), 
                       include.global = T, symmetry.test = T, normality.test = F)
df <- plyr::rbind.fill(
  df, do.call(plyr::rbind.fill, lapply(lfatores2, FUN = function(f) {
    if (nrow(dat) > 0 && sum(!is.na(unique(dat[[f]]))) > 1)
      get.descriptives(dat, c(dv.pre,dv.pos), c("grupo", f),
                       symmetry.test = T, normality.test = F)
    }))
)
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning:
## ! There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning in `stats::qt()`:
## ! NaNs produced
## There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning:
## ! There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning in `stats::qt()`:
## ! NaNs produced
df <- df[,c(fatores1[fatores1 %in% colnames(df)],"variable",
            colnames(df)[!colnames(df) %in% c(fatores1,"variable")])]
grupo genero zona.participante zona.escola score.vocab.ensinado.quintile variable n mean median min max sd se ci iqr symmetry skewness kurtosis
Controle score.vocab.ensinado.pre 50 2.420 2.0 0 6 1.527 0.216 0.434 2.75 YES 0.267 -0.595
Experimental score.vocab.ensinado.pre 42 2.833 3.0 0 7 1.447 0.223 0.451 2.00 YES 0.378 -0.084
score.vocab.ensinado.pre 92 2.609 3.0 0 7 1.497 0.156 0.310 3.00 YES 0.288 -0.299
Controle score.vocab.ensinado.pos 50 2.360 2.0 0 8 1.711 0.242 0.486 2.00 NO 1.048 1.269
Experimental score.vocab.ensinado.pos 42 2.643 2.0 0 6 1.694 0.261 0.528 2.75 YES 0.142 -0.982
score.vocab.ensinado.pos 92 2.489 2.0 0 8 1.700 0.177 0.352 3.00 NO 0.642 0.178
Controle F score.vocab.ensinado.pre 25 2.640 3.0 0 6 1.497 0.299 0.618 2.00 YES 0.098 -0.588
Controle M score.vocab.ensinado.pre 25 2.200 2.0 0 6 1.555 0.311 0.642 2.00 YES 0.447 -0.620
Experimental F score.vocab.ensinado.pre 19 2.579 3.0 0 5 1.427 0.327 0.688 2.00 YES 0.067 -1.086
Experimental M score.vocab.ensinado.pre 23 3.043 3.0 1 7 1.461 0.305 0.632 2.00 NO 0.598 0.121
Controle F score.vocab.ensinado.pos 25 2.400 2.0 0 8 1.871 0.374 0.772 2.00 NO 1.415 1.611
Controle M score.vocab.ensinado.pos 25 2.320 2.0 0 6 1.574 0.315 0.650 2.00 YES 0.289 -0.488
Experimental F score.vocab.ensinado.pos 19 2.579 2.0 0 6 1.835 0.421 0.885 3.00 YES 0.136 -1.266
Experimental M score.vocab.ensinado.pos 23 2.696 2.0 0 6 1.608 0.335 0.695 2.00 YES 0.163 -0.904
Controle Rural score.vocab.ensinado.pre 12 2.000 2.5 0 4 1.537 0.444 0.977 2.25 YES -0.138 -1.717
Controle Urbana score.vocab.ensinado.pre 26 2.500 2.0 0 6 1.631 0.320 0.659 3.00 YES 0.479 -0.610
Controle score.vocab.ensinado.pre 12 2.667 3.0 1 5 1.303 0.376 0.828 1.50 YES 0.117 -1.277
Experimental Rural score.vocab.ensinado.pre 15 2.800 3.0 1 5 1.265 0.327 0.700 1.50 YES 0.348 -1.072
Experimental Urbana score.vocab.ensinado.pre 16 2.938 3.0 0 7 1.611 0.403 0.859 2.00 NO 0.543 0.443
Experimental score.vocab.ensinado.pre 11 2.727 3.0 1 5 1.555 0.469 1.045 3.00 YES -0.024 -1.864
Controle Rural score.vocab.ensinado.pos 12 2.833 2.0 1 8 1.801 0.520 1.144 1.00 NO 1.854 2.776
Controle Urbana score.vocab.ensinado.pos 26 2.192 2.0 0 6 1.877 0.368 0.758 2.00 NO 0.780 -0.425
Controle score.vocab.ensinado.pos 12 2.250 2.0 0 4 1.215 0.351 0.772 1.25 YES -0.157 -1.082
Experimental Rural score.vocab.ensinado.pos 15 2.333 2.0 0 5 1.496 0.386 0.828 2.00 YES 0.062 -1.140
Experimental Urbana score.vocab.ensinado.pos 16 3.000 3.5 0 6 1.826 0.456 0.973 2.25 YES -0.246 -1.234
Experimental score.vocab.ensinado.pos 11 2.545 2.0 0 6 1.809 0.545 1.215 2.00 NO 0.532 -0.989
Controle Rural score.vocab.ensinado.pre 14 2.214 2.0 0 5 1.528 0.408 0.882 2.00 YES 0.267 -1.439
Controle Urbana score.vocab.ensinado.pre 36 2.500 2.0 0 6 1.540 0.257 0.521 2.25 YES 0.251 -0.415
Experimental Rural score.vocab.ensinado.pre 13 2.846 3.0 1 5 1.144 0.317 0.691 2.00 YES 0.272 -1.155
Experimental Urbana score.vocab.ensinado.pre 29 2.828 3.0 0 7 1.583 0.294 0.602 2.00 YES 0.378 -0.272
Controle Rural score.vocab.ensinado.pos 14 2.929 2.5 1 8 1.900 0.508 1.097 1.00 NO 1.469 1.286
Controle Urbana score.vocab.ensinado.pos 36 2.139 2.0 0 6 1.606 0.268 0.543 2.00 NO 0.665 -0.090
Experimental Rural score.vocab.ensinado.pos 13 2.538 2.0 0 6 1.941 0.538 1.173 3.00 YES 0.356 -1.288
Experimental Urbana score.vocab.ensinado.pos 29 2.690 3.0 0 6 1.606 0.298 0.611 2.00 YES -0.008 -0.968
Controle 1st quintile score.vocab.ensinado.pre 16 0.688 1.0 0 1 0.479 0.120 0.255 1.00 few data 0.000 0.000
Controle 2nd quintile score.vocab.ensinado.pre 21 2.524 3.0 2 3 0.512 0.112 0.233 1.00 few data 0.000 0.000
Controle 3rd quintile score.vocab.ensinado.pre 13 4.385 4.0 4 6 0.768 0.213 0.464 0.00 few data 0.000 0.000
Experimental 1st quintile score.vocab.ensinado.pre 8 0.875 1.0 0 1 0.354 0.125 0.296 0.00 few data 0.000 0.000
Experimental 2nd quintile score.vocab.ensinado.pre 19 2.421 2.0 2 3 0.507 0.116 0.244 1.00 few data 0.000 0.000
Experimental 3rd quintile score.vocab.ensinado.pre 14 4.214 4.0 4 5 0.426 0.114 0.246 0.00 few data 0.000 0.000
Experimental 4th quintile score.vocab.ensinado.pre 1 7.000 7.0 7 7 0.00 few data 0.000 0.000
Controle 1st quintile score.vocab.ensinado.pos 16 2.000 2.0 0 6 1.549 0.387 0.826 1.25 NO 0.807 0.516
Controle 2nd quintile score.vocab.ensinado.pos 21 2.429 2.0 0 8 1.859 0.406 0.846 2.00 NO 1.352 1.858
Controle 3rd quintile score.vocab.ensinado.pos 13 2.692 2.0 0 6 1.702 0.472 1.029 2.00 YES 0.351 -0.949
Experimental 1st quintile score.vocab.ensinado.pos 8 2.750 2.0 1 5 1.389 0.491 1.161 2.00 YES 0.385 -1.617
Experimental 2nd quintile score.vocab.ensinado.pos 19 3.158 3.0 0 6 1.642 0.377 0.791 2.00 YES -0.025 -0.863
Experimental 3rd quintile score.vocab.ensinado.pos 14 1.786 1.5 0 5 1.718 0.459 0.992 2.50 NO 0.645 -0.870
Experimental 4th quintile score.vocab.ensinado.pos 1 4.000 4.0 4 4 0.00 few data 0.000 0.000

ANCOVA and Pairwise for one factor: grupo

Without remove non-normal data

pdat = remove_group_data(dat[!is.na(dat[["grupo"]]),], "score.vocab.ensinado.pos", "grupo")

pdat.long <- rbind(pdat[,c("id","grupo")], pdat[,c("id","grupo")])
pdat.long[["time"]] <- c(rep("pre", nrow(pdat)), rep("pos", nrow(pdat)))
pdat.long[["time"]] <- factor(pdat.long[["time"]], c("pre","pos"))
pdat.long[["score.vocab.ensinado"]] <- c(pdat[["score.vocab.ensinado.pre"]], pdat[["score.vocab.ensinado.pos"]])

aov = anova_test(pdat, score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo)
laov[["grupo"]] <- get_anova_table(aov)
pwc <- emmeans_test(pdat, score.vocab.ensinado.pos ~ grupo, covariate = score.vocab.ensinado.pre,
                    p.adjust.method = "bonferroni")
pwc.long <- emmeans_test(dplyr::group_by_at(pdat.long, "grupo"),
                          score.vocab.ensinado ~ time,
                          p.adjust.method = "bonferroni")
lpwc[["grupo"]] <- plyr::rbind.fill(pwc, pwc.long)
ds <- get.descriptives(pdat, "score.vocab.ensinado.pos", "grupo", covar = "score.vocab.ensinado.pre")
ds <- merge(ds[ds$variable != "score.vocab.ensinado.pre",],
            ds[ds$variable == "score.vocab.ensinado.pre", !colnames(ds) %in% c("variable")],
            by = "grupo", all.x = T, suffixes = c("", ".score.vocab.ensinado.pre"))
ds <- merge(get_emmeans(pwc), ds, by = "grupo", suffixes = c(".emms", ""))
ds <- ds[,c("grupo","n","mean.score.vocab.ensinado.pre","se.score.vocab.ensinado.pre","mean","se",
            "emmean","se.emms","conf.low","conf.high")]

colnames(ds) <- c("grupo", "N", paste0(c("M","SE")," (pre)"),
                  paste0(c("M","SE"), " (unadj)"),
                  paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")

lemms[["grupo"]] <- ds

Computing ANCOVA and PairWise After removing non-normal data (OK)

wdat = pdat 

res = residuals(lm(score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo, data = wdat))
non.normal = getNonNormal(res, wdat$id, plimit = 0.05)

wdat = wdat[!wdat$id %in% non.normal,]

wdat.long <- rbind(wdat[,c("id","grupo")], wdat[,c("id","grupo")])
wdat.long[["time"]] <- c(rep("pre", nrow(wdat)), rep("pos", nrow(wdat)))
wdat.long[["time"]] <- factor(wdat.long[["time"]], c("pre","pos"))
wdat.long[["score.vocab.ensinado"]] <- c(wdat[["score.vocab.ensinado.pre"]], wdat[["score.vocab.ensinado.pos"]])

ldat[["grupo"]] = wdat

(non.normal)
## [1] "P264" "P286" "P78"  "P97"  "P272" "P85"  "P273" "P82"
aov = anova_test(wdat, score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo)
laov[["grupo"]] <- merge(get_anova_table(aov), laov[["grupo"]],
                            by="Effect", suffixes = c("","'"))

(df = get_anova_table(aov))
## ANOVA Table (type II tests)
## 
##                     Effect DFn DFd     F     p p<.05   ges
## 1 score.vocab.ensinado.pre   1  81 0.429 0.514       0.005
## 2                    grupo   1  81 3.174 0.079       0.038
Effect DFn DFd F p p<.05 ges
score.vocab.ensinado.pre 1 81 0.429 0.514 0.005
grupo 1 81 3.174 0.079 0.038
pwc <- emmeans_test(wdat, score.vocab.ensinado.pos ~ grupo, covariate = score.vocab.ensinado.pre,
                    p.adjust.method = "bonferroni")
term .y. group1 group2 df statistic p p.adj p.adj.signif
score.vocab.ensinado.pre*grupo score.vocab.ensinado.pos Controle Experimental 81 -1.782 0.079 0.079 ns
pwc.long <- emmeans_test(dplyr::group_by_at(wdat.long, "grupo"),
                         score.vocab.ensinado ~ time,
                         p.adjust.method = "bonferroni")
lpwc[["grupo"]] <- merge(plyr::rbind.fill(pwc, pwc.long), lpwc[["grupo"]],
                            by=c("grupo","term",".y.","group1","group2"),
                            suffixes = c("","'"))
grupo term .y. group1 group2 df statistic p p.adj p.adj.signif
Controle time score.vocab.ensinado pre pos 164 0.567 0.571 0.571 ns
Experimental time score.vocab.ensinado pre pos 164 -0.246 0.806 0.806 ns
ds <- get.descriptives(wdat, "score.vocab.ensinado.pos", "grupo", covar = "score.vocab.ensinado.pre")
ds <- merge(ds[ds$variable != "score.vocab.ensinado.pre",],
            ds[ds$variable == "score.vocab.ensinado.pre", !colnames(ds) %in% c("variable")],
            by = "grupo", all.x = T, suffixes = c("", ".score.vocab.ensinado.pre"))
ds <- merge(get_emmeans(pwc), ds, by = "grupo", suffixes = c(".emms", ""))
ds <- ds[,c("grupo","n","mean.score.vocab.ensinado.pre","se.score.vocab.ensinado.pre","mean","se",
            "emmean","se.emms","conf.low","conf.high")]

colnames(ds) <- c("grupo", "N", paste0(c("M","SE")," (pre)"),
                  paste0(c("M","SE"), " (unadj)"),
                  paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")

lemms[["grupo"]] <- merge(ds, lemms[["grupo"]], by=c("grupo"), suffixes = c("","'"))
grupo N M (pre) SE (pre) M (unadj) SE (unadj) M (adj) SE (adj) conf.low conf.high
Controle 48 2.333 0.211 2.167 0.207 2.177 0.205 1.769 2.585
Experimental 36 2.667 0.246 2.750 0.230 2.737 0.237 2.265 3.208

Plots for ancova

plots <- oneWayAncovaPlots(
  wdat, "score.vocab.ensinado.pos", "grupo", aov, list("grupo"=pwc), addParam = c("mean_ci"),
  font.label.size=10, step.increase=0.05, p.label="p.adj",
  subtitle = which(aov$Effect == "grupo"))
if (!is.null(nrow(plots[["grupo"]]$data)))
  plots[["grupo"]] + ggplot2::ylab("Vocabulary taught") + 
  if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)

plots <- oneWayAncovaBoxPlots(
  wdat, "score.vocab.ensinado.pos", "grupo", aov, pwc, covar = "score.vocab.ensinado.pre",
  theme = "classic", color = color[["grupo"]],
  subtitle = which(aov$Effect == "grupo"))
if (length(unique(wdat[["grupo"]])) > 1)
  plots[["grupo"]] + ggplot2::ylab("Vocabulary taught") +
  ggplot2::scale_x_discrete(labels=c('pre', 'pos')) +
  if (ymin < ymax) ggplot2::ylim(ymin, ymax)

if (length(unique(wdat.long[["grupo"]])) > 1)
  plots <- oneWayAncovaBoxPlots(
    wdat.long, "score.vocab.ensinado", "grupo", aov, pwc.long,
    pre.post = "time", theme = "classic", color = color$prepost)
if (length(unique(wdat.long[["grupo"]])) > 1)
  plots[["grupo"]] + ggplot2::ylab("Vocabulary taught") +
  if (ymin < ymax) ggplot2::ylim(ymin, ymax) 

Checking linearity assumption

ggscatter(wdat, x = "score.vocab.ensinado.pre", y = "score.vocab.ensinado.pos", size = 0.5,
          color = "grupo", add = "reg.line")+
  stat_regline_equation(
    aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = grupo)
  ) +
  ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo"))) +
  ggplot2::scale_color_manual(values = color[["grupo"]]) +
  ggplot2::ylab("Vocabulary taught")  +
  if (ymin < ymax) ggplot2::ylim(ymin, ymax)

Checking normality and homogeneity

res <- augment(lm(score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo, data = wdat))
shapiro_test(res$.resid)
## # A tibble: 1 × 3
##   variable   statistic p.value
##   <chr>          <dbl>   <dbl>
## 1 res$.resid     0.977   0.145
levene_test(res, .resid ~ grupo)
## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     1    82     0.757 0.387

ANCOVA and Pairwise for two factors grupo:genero

Without remove non-normal data

pdat = remove_group_data(dat[!is.na(dat[["grupo"]]) & !is.na(dat[["genero"]]),],
                         "score.vocab.ensinado.pos", c("grupo","genero"))
pdat = pdat[pdat[["genero"]] %in% do.call(
  intersect, lapply(unique(pdat[["grupo"]]), FUN = function(x) {
    unique(pdat[["genero"]][which(pdat[["grupo"]] == x)])
  })),]
pdat[["grupo"]] = factor(pdat[["grupo"]], level[["grupo"]])
pdat[["genero"]] = factor(
  pdat[["genero"]],
  level[["genero"]][level[["genero"]] %in% unique(pdat[["genero"]])])

pdat.long <- rbind(pdat[,c("id","grupo","genero")], pdat[,c("id","grupo","genero")])
pdat.long[["time"]] <- c(rep("pre", nrow(pdat)), rep("pos", nrow(pdat)))
pdat.long[["time"]] <- factor(pdat.long[["time"]], c("pre","pos"))
pdat.long[["score.vocab.ensinado"]] <- c(pdat[["score.vocab.ensinado.pre"]], pdat[["score.vocab.ensinado.pos"]])

if (length(unique(pdat[["genero"]])) >= 2) {
  aov = anova_test(pdat, score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo*genero)
  laov[["grupo:genero"]] <- get_anova_table(aov)
}
if (length(unique(pdat[["genero"]])) >= 2) {
  pwcs <- list()
  pwcs[["genero"]] <- emmeans_test(
    group_by(pdat, grupo), score.vocab.ensinado.pos ~ genero,
    covariate = score.vocab.ensinado.pre, p.adjust.method = "bonferroni")
  pwcs[["grupo"]] <- emmeans_test(
    group_by(pdat, genero), score.vocab.ensinado.pos ~ grupo,
    covariate = score.vocab.ensinado.pre, p.adjust.method = "bonferroni")
  
  pwc <- plyr::rbind.fill(pwcs[["grupo"]], pwcs[["genero"]])
  pwc <- pwc[,c("grupo","genero", colnames(pwc)[!colnames(pwc) %in% c("grupo","genero")])]
}
if (length(unique(pdat[["genero"]])) >= 2) {
  pwc.long <- emmeans_test(dplyr::group_by_at(pdat.long, c("grupo","genero")),
                           score.vocab.ensinado ~ time,
                           p.adjust.method = "bonferroni")
  lpwc[["grupo:genero"]] <- plyr::rbind.fill(pwc, pwc.long)
}
if (length(unique(pdat[["genero"]])) >= 2) {
  ds <- get.descriptives(pdat, "score.vocab.ensinado.pos", c("grupo","genero"), covar = "score.vocab.ensinado.pre")
  ds <- merge(ds[ds$variable != "score.vocab.ensinado.pre",],
              ds[ds$variable == "score.vocab.ensinado.pre", !colnames(ds) %in% c("variable")],
              by = c("grupo","genero"), all.x = T, suffixes = c("", ".score.vocab.ensinado.pre"))
  ds <- merge(get_emmeans(pwcs[["grupo"]]), ds,
              by = c("grupo","genero"), suffixes = c(".emms", ""))
  ds <- ds[,c("grupo","genero","n","mean.score.vocab.ensinado.pre","se.score.vocab.ensinado.pre","mean","se",
              "emmean","se.emms","conf.low","conf.high")]
  
  colnames(ds) <- c("grupo","genero", "N", paste0(c("M","SE")," (pre)"),
                    paste0(c("M","SE"), " (unadj)"),
                    paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")
  
  lemms[["grupo:genero"]] <- ds
}

Computing ANCOVA and PairWise After removing non-normal data (OK)

if (length(unique(pdat[["genero"]])) >= 2) {
  wdat = pdat 
  
  res = residuals(lm(score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo*genero, data = wdat))
  non.normal = getNonNormal(res, wdat$id, plimit = 0.05)
  
  wdat = wdat[!wdat$id %in% non.normal,]
  
  wdat.long <- rbind(wdat[,c("id","grupo","genero")], wdat[,c("id","grupo","genero")])
  wdat.long[["time"]] <- c(rep("pre", nrow(wdat)), rep("pos", nrow(wdat)))
  wdat.long[["time"]] <- factor(wdat.long[["time"]], c("pre","pos"))
  wdat.long[["score.vocab.ensinado"]] <- c(wdat[["score.vocab.ensinado.pre"]], wdat[["score.vocab.ensinado.pos"]])
  
  
  ldat[["grupo:genero"]] = wdat
  
  (non.normal)
}
## [1] "P97"  "P286" "P272" "P264" "P78"  "P85"  "P273" "P82"
if (length(unique(pdat[["genero"]])) >= 2) {
  aov = anova_test(wdat, score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo*genero)
  laov[["grupo:genero"]] <- merge(get_anova_table(aov), laov[["grupo:genero"]],
                                         by="Effect", suffixes = c("","'"))
  df = get_anova_table(aov)
}
Effect DFn DFd F p p<.05 ges
score.vocab.ensinado.pre 1 79 0.394 0.532 0.005
grupo 1 79 3.053 0.084 0.037
genero 1 79 0.028 0.868 0.000
grupo:genero 1 79 0.002 0.966 0.000
if (length(unique(pdat[["genero"]])) >= 2) {
  pwcs <- list()
  pwcs[["genero"]] <- emmeans_test(
    group_by(wdat, grupo), score.vocab.ensinado.pos ~ genero,
    covariate = score.vocab.ensinado.pre, p.adjust.method = "bonferroni")
  pwcs[["grupo"]] <- emmeans_test(
    group_by(wdat, genero), score.vocab.ensinado.pos ~ grupo,
    covariate = score.vocab.ensinado.pre, p.adjust.method = "bonferroni")
  
  pwc <- plyr::rbind.fill(pwcs[["grupo"]], pwcs[["genero"]])
  pwc <- pwc[,c("grupo","genero", colnames(pwc)[!colnames(pwc) %in% c("grupo","genero")])]
}
grupo genero term .y. group1 group2 df statistic p p.adj p.adj.signif
F score.vocab.ensinado.pre*grupo score.vocab.ensinado.pos Controle Experimental 79 -1.169 0.246 0.246 ns
M score.vocab.ensinado.pre*grupo score.vocab.ensinado.pos Controle Experimental 79 -1.280 0.204 0.204 ns
Controle score.vocab.ensinado.pre*genero score.vocab.ensinado.pos F M 79 -0.096 0.923 0.923 ns
Experimental score.vocab.ensinado.pre*genero score.vocab.ensinado.pos F M 79 -0.141 0.889 0.889 ns
if (length(unique(pdat[["genero"]])) >= 2) {
  pwc.long <- emmeans_test(dplyr::group_by_at(wdat.long, c("grupo","genero")),
                           score.vocab.ensinado ~ time,
                           p.adjust.method = "bonferroni")
  lpwc[["grupo:genero"]] <- merge(plyr::rbind.fill(pwc, pwc.long),
                                         lpwc[["grupo:genero"]],
                                         by=c("grupo","genero","term",".y.","group1","group2"),
                                         suffixes = c("","'"))
}
grupo genero term .y. group1 group2 df statistic p p.adj p.adj.signif
Controle F time score.vocab.ensinado pre pos 160 1.102 0.272 0.272 ns
Controle M time score.vocab.ensinado pre pos 160 -0.300 0.764 0.764 ns
Experimental F time score.vocab.ensinado pre pos 160 -0.736 0.463 0.463 ns
Experimental M time score.vocab.ensinado pre pos 160 0.329 0.742 0.742 ns
if (length(unique(pdat[["genero"]])) >= 2) {
  ds <- get.descriptives(wdat, "score.vocab.ensinado.pos", c("grupo","genero"), covar = "score.vocab.ensinado.pre")
  ds <- merge(ds[ds$variable != "score.vocab.ensinado.pre",],
              ds[ds$variable == "score.vocab.ensinado.pre", !colnames(ds) %in% c("variable")],
              by = c("grupo","genero"), all.x = T, suffixes = c("", ".score.vocab.ensinado.pre"))
  ds <- merge(get_emmeans(pwcs[["grupo"]]), ds,
              by = c("grupo","genero"), suffixes = c(".emms", ""))
  ds <- ds[,c("grupo","genero","n","mean.score.vocab.ensinado.pre","se.score.vocab.ensinado.pre",
              "mean","se","emmean","se.emms","conf.low","conf.high")]
  
  colnames(ds) <- c("grupo","genero", "N", paste0(c("M","SE")," (pre)"),
                    paste0(c("M","SE"), " (unadj)"),
                    paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")
  
  lemms[["grupo:genero"]] <- merge(ds, lemms[["grupo:genero"]],
                                          by=c("grupo","genero"), suffixes = c("","'"))
}
grupo genero N M (pre) SE (pre) M (unadj) SE (unadj) M (adj) SE (adj) conf.low conf.high
Controle F 24 2.625 0.312 2.167 0.305 2.156 0.293 1.573 2.740
Controle M 24 2.042 0.279 2.167 0.287 2.197 0.297 1.606 2.787
Experimental F 16 2.312 0.338 2.688 0.384 2.699 0.359 1.984 3.413
Experimental M 20 2.950 0.344 2.800 0.287 2.767 0.325 2.120 3.414

Plots for ancova

if (length(unique(pdat[["genero"]])) >= 2) {
  ggPlotAoC2(pwcs, "grupo", "genero", aov, ylab = "Vocabulary taught",
             subtitle = which(aov$Effect == "grupo:genero"), addParam = "errorbar") +
    ggplot2::scale_color_manual(values = color[["genero"]]) +
    ggplot2::ylab("Vocabulary taught") +
    if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

if (length(unique(pdat[["genero"]])) >= 2) {
  ggPlotAoC2(pwcs, "genero", "grupo", aov, ylab = "Vocabulary taught",
               subtitle = which(aov$Effect == "grupo:genero"), addParam = "errorbar") +
    ggplot2::scale_color_manual(values = color[["grupo"]]) +
    ggplot2::ylab("Vocabulary taught") +
    if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

if (length(unique(pdat[["genero"]])) >= 2) {
  plots <- twoWayAncovaBoxPlots(
    wdat, "score.vocab.ensinado.pos", c("grupo","genero"), aov, pwcs, covar = "score.vocab.ensinado.pre",
    theme = "classic", color = color[["grupo:genero"]],
    subtitle = which(aov$Effect == "grupo:genero"))
}
if (length(unique(pdat[["genero"]])) >= 2) {
  plots[["grupo:genero"]] + ggplot2::ylab("Vocabulary taught") +
  ggplot2::scale_x_discrete(labels=c('pre', 'pos')) +
  if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}
## Warning: No shared levels found between `names(values)` of the manual scale and the data's colour values.

if (length(unique(pdat[["genero"]])) >= 2) {
  plots <- twoWayAncovaBoxPlots(
    wdat.long, "score.vocab.ensinado", c("grupo","genero"), aov, pwc.long,
    pre.post = "time",
    theme = "classic", color = color$prepost)
}
if (length(unique(pdat[["genero"]])) >= 2) 
  plots[["grupo:genero"]] + ggplot2::ylab("Vocabulary taught") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)

Checking linearity assumption

if (length(unique(pdat[["genero"]])) >= 2) {
  ggscatter(wdat, x = "score.vocab.ensinado.pre", y = "score.vocab.ensinado.pos", size = 0.5,
            facet.by = c("grupo","genero"), add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"))
    ) + ggplot2::ylab("Vocabulary taught") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["genero"]])) >= 2) {
  ggscatter(wdat, x = "score.vocab.ensinado.pre", y = "score.vocab.ensinado.pos", size = 0.5,
            color = "grupo", facet.by = "genero", add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = grupo)
    ) +
    ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo:genero"))) +
    ggplot2::scale_color_manual(values = color[["grupo"]]) +
    ggplot2::ylab("Vocabulary taught") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["genero"]])) >= 2) {
  ggscatter(wdat, x = "score.vocab.ensinado.pre", y = "score.vocab.ensinado.pos", size = 0.5,
            color = "genero", facet.by = "grupo", add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = genero)
    ) +
    ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo:genero"))) +
    ggplot2::scale_color_manual(values = color[["genero"]]) +
    ggplot2::ylab("Vocabulary taught") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

Checking normality and homogeneity

if (length(unique(pdat[["genero"]])) >= 2) 
  res <- augment(lm(score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo*genero, data = wdat))
if (length(unique(pdat[["genero"]])) >= 2)
  shapiro_test(res$.resid)
## # A tibble: 1 × 3
##   variable   statistic p.value
##   <chr>          <dbl>   <dbl>
## 1 res$.resid     0.977   0.133
if (length(unique(pdat[["genero"]])) >= 2) 
  levene_test(res, .resid ~ grupo*genero)
## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     3    80     0.439 0.726

ANCOVA and Pairwise for two factors grupo:zona.participante

Without remove non-normal data

pdat = remove_group_data(dat[!is.na(dat[["grupo"]]) & !is.na(dat[["zona.participante"]]),],
                         "score.vocab.ensinado.pos", c("grupo","zona.participante"))
pdat = pdat[pdat[["zona.participante"]] %in% do.call(
  intersect, lapply(unique(pdat[["grupo"]]), FUN = function(x) {
    unique(pdat[["zona.participante"]][which(pdat[["grupo"]] == x)])
  })),]
pdat[["grupo"]] = factor(pdat[["grupo"]], level[["grupo"]])
pdat[["zona.participante"]] = factor(
  pdat[["zona.participante"]],
  level[["zona.participante"]][level[["zona.participante"]] %in% unique(pdat[["zona.participante"]])])

pdat.long <- rbind(pdat[,c("id","grupo","zona.participante")], pdat[,c("id","grupo","zona.participante")])
pdat.long[["time"]] <- c(rep("pre", nrow(pdat)), rep("pos", nrow(pdat)))
pdat.long[["time"]] <- factor(pdat.long[["time"]], c("pre","pos"))
pdat.long[["score.vocab.ensinado"]] <- c(pdat[["score.vocab.ensinado.pre"]], pdat[["score.vocab.ensinado.pos"]])

if (length(unique(pdat[["zona.participante"]])) >= 2) {
  aov = anova_test(pdat, score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo*zona.participante)
  laov[["grupo:zona.participante"]] <- get_anova_table(aov)
}
if (length(unique(pdat[["zona.participante"]])) >= 2) {
  pwcs <- list()
  pwcs[["zona.participante"]] <- emmeans_test(
    group_by(pdat, grupo), score.vocab.ensinado.pos ~ zona.participante,
    covariate = score.vocab.ensinado.pre, p.adjust.method = "bonferroni")
  pwcs[["grupo"]] <- emmeans_test(
    group_by(pdat, zona.participante), score.vocab.ensinado.pos ~ grupo,
    covariate = score.vocab.ensinado.pre, p.adjust.method = "bonferroni")
  
  pwc <- plyr::rbind.fill(pwcs[["grupo"]], pwcs[["zona.participante"]])
  pwc <- pwc[,c("grupo","zona.participante", colnames(pwc)[!colnames(pwc) %in% c("grupo","zona.participante")])]
}
if (length(unique(pdat[["zona.participante"]])) >= 2) {
  pwc.long <- emmeans_test(dplyr::group_by_at(pdat.long, c("grupo","zona.participante")),
                           score.vocab.ensinado ~ time,
                           p.adjust.method = "bonferroni")
  lpwc[["grupo:zona.participante"]] <- plyr::rbind.fill(pwc, pwc.long)
}
if (length(unique(pdat[["zona.participante"]])) >= 2) {
  ds <- get.descriptives(pdat, "score.vocab.ensinado.pos", c("grupo","zona.participante"), covar = "score.vocab.ensinado.pre")
  ds <- merge(ds[ds$variable != "score.vocab.ensinado.pre",],
              ds[ds$variable == "score.vocab.ensinado.pre", !colnames(ds) %in% c("variable")],
              by = c("grupo","zona.participante"), all.x = T, suffixes = c("", ".score.vocab.ensinado.pre"))
  ds <- merge(get_emmeans(pwcs[["grupo"]]), ds,
              by = c("grupo","zona.participante"), suffixes = c(".emms", ""))
  ds <- ds[,c("grupo","zona.participante","n","mean.score.vocab.ensinado.pre","se.score.vocab.ensinado.pre","mean","se",
              "emmean","se.emms","conf.low","conf.high")]
  
  colnames(ds) <- c("grupo","zona.participante", "N", paste0(c("M","SE")," (pre)"),
                    paste0(c("M","SE"), " (unadj)"),
                    paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")
  
  lemms[["grupo:zona.participante"]] <- ds
}

Computing ANCOVA and PairWise After removing non-normal data (OK)

if (length(unique(pdat[["zona.participante"]])) >= 2) {
  wdat = pdat 
  
  res = residuals(lm(score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo*zona.participante, data = wdat))
  non.normal = getNonNormal(res, wdat$id, plimit = 0.05)
  
  wdat = wdat[!wdat$id %in% non.normal,]
  
  wdat.long <- rbind(wdat[,c("id","grupo","zona.participante")], wdat[,c("id","grupo","zona.participante")])
  wdat.long[["time"]] <- c(rep("pre", nrow(wdat)), rep("pos", nrow(wdat)))
  wdat.long[["time"]] <- factor(wdat.long[["time"]], c("pre","pos"))
  wdat.long[["score.vocab.ensinado"]] <- c(wdat[["score.vocab.ensinado.pre"]], wdat[["score.vocab.ensinado.pos"]])
  
  
  ldat[["grupo:zona.participante"]] = wdat
  
  (non.normal)
}
## NULL
if (length(unique(pdat[["zona.participante"]])) >= 2) {
  aov = anova_test(wdat, score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo*zona.participante)
  laov[["grupo:zona.participante"]] <- merge(get_anova_table(aov), laov[["grupo:zona.participante"]],
                                         by="Effect", suffixes = c("","'"))
  df = get_anova_table(aov)
}
Effect DFn DFd F p p<.05 ges
score.vocab.ensinado.pre 1 64 1.945 0.168 0.029
grupo 1 64 0.155 0.695 0.002
zona.participante 1 64 0.022 0.881 0.000
grupo:zona.participante 1 64 2.431 0.124 0.037
if (length(unique(pdat[["zona.participante"]])) >= 2) {
  pwcs <- list()
  pwcs[["zona.participante"]] <- emmeans_test(
    group_by(wdat, grupo), score.vocab.ensinado.pos ~ zona.participante,
    covariate = score.vocab.ensinado.pre, p.adjust.method = "bonferroni")
  pwcs[["grupo"]] <- emmeans_test(
    group_by(wdat, zona.participante), score.vocab.ensinado.pos ~ grupo,
    covariate = score.vocab.ensinado.pre, p.adjust.method = "bonferroni")
  
  pwc <- plyr::rbind.fill(pwcs[["grupo"]], pwcs[["zona.participante"]])
  pwc <- pwc[,c("grupo","zona.participante", colnames(pwc)[!colnames(pwc) %in% c("grupo","zona.participante")])]
}
grupo zona.participante term .y. group1 group2 df statistic p p.adj p.adj.signif
Rural score.vocab.ensinado.pre*grupo score.vocab.ensinado.pos Controle Experimental 64 0.951 0.345 0.345 ns
Urbana score.vocab.ensinado.pre*grupo score.vocab.ensinado.pos Controle Experimental 64 -1.279 0.206 0.206 ns
Controle score.vocab.ensinado.pre*zona.participante score.vocab.ensinado.pos Rural Urbana 64 1.195 0.237 0.237 ns
Experimental score.vocab.ensinado.pre*zona.participante score.vocab.ensinado.pos Rural Urbana 64 -1.008 0.317 0.317 ns
if (length(unique(pdat[["zona.participante"]])) >= 2) {
  pwc.long <- emmeans_test(dplyr::group_by_at(wdat.long, c("grupo","zona.participante")),
                           score.vocab.ensinado ~ time,
                           p.adjust.method = "bonferroni")
  lpwc[["grupo:zona.participante"]] <- merge(plyr::rbind.fill(pwc, pwc.long),
                                         lpwc[["grupo:zona.participante"]],
                                         by=c("grupo","zona.participante","term",".y.","group1","group2"),
                                         suffixes = c("","'"))
}
grupo zona.participante term .y. group1 group2 df statistic p p.adj p.adj.signif
Controle Rural time score.vocab.ensinado pre pos 130 -1.229 0.221 0.221 ns
Controle Urbana time score.vocab.ensinado pre pos 130 0.668 0.506 0.506 ns
Experimental Rural time score.vocab.ensinado pre pos 130 0.769 0.443 0.443 ns
Experimental Urbana time score.vocab.ensinado pre pos 130 -0.106 0.915 0.915 ns
if (length(unique(pdat[["zona.participante"]])) >= 2) {
  ds <- get.descriptives(wdat, "score.vocab.ensinado.pos", c("grupo","zona.participante"), covar = "score.vocab.ensinado.pre")
  ds <- merge(ds[ds$variable != "score.vocab.ensinado.pre",],
              ds[ds$variable == "score.vocab.ensinado.pre", !colnames(ds) %in% c("variable")],
              by = c("grupo","zona.participante"), all.x = T, suffixes = c("", ".score.vocab.ensinado.pre"))
  ds <- merge(get_emmeans(pwcs[["grupo"]]), ds,
              by = c("grupo","zona.participante"), suffixes = c(".emms", ""))
  ds <- ds[,c("grupo","zona.participante","n","mean.score.vocab.ensinado.pre","se.score.vocab.ensinado.pre",
              "mean","se","emmean","se.emms","conf.low","conf.high")]
  
  colnames(ds) <- c("grupo","zona.participante", "N", paste0(c("M","SE")," (pre)"),
                    paste0(c("M","SE"), " (unadj)"),
                    paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")
  
  lemms[["grupo:zona.participante"]] <- merge(ds, lemms[["grupo:zona.participante"]],
                                          by=c("grupo","zona.participante"), suffixes = c("","'"))
}
grupo zona.participante N M (pre) SE (pre) M (unadj) SE (unadj) M (adj) SE (adj) conf.low conf.high
Controle Rural 12 2.000 0.444 2.833 0.520 2.948 0.516 1.918 3.978
Controle Urbana 26 2.500 0.320 2.192 0.368 2.208 0.346 1.517 2.899
Experimental Rural 15 2.800 0.327 2.333 0.386 2.290 0.456 1.378 3.201
Experimental Urbana 16 2.938 0.403 3.000 0.456 2.929 0.444 2.043 3.816

Plots for ancova

if (length(unique(pdat[["zona.participante"]])) >= 2) {
  ggPlotAoC2(pwcs, "grupo", "zona.participante", aov, ylab = "Vocabulary taught",
             subtitle = which(aov$Effect == "grupo:zona.participante"), addParam = "errorbar") +
    ggplot2::scale_color_manual(values = color[["zona.participante"]]) +
    ggplot2::ylab("Vocabulary taught") +
    if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

if (length(unique(pdat[["zona.participante"]])) >= 2) {
  ggPlotAoC2(pwcs, "zona.participante", "grupo", aov, ylab = "Vocabulary taught",
               subtitle = which(aov$Effect == "grupo:zona.participante"), addParam = "errorbar") +
    ggplot2::scale_color_manual(values = color[["grupo"]]) +
    ggplot2::ylab("Vocabulary taught") +
    if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

if (length(unique(pdat[["zona.participante"]])) >= 2) {
  plots <- twoWayAncovaBoxPlots(
    wdat, "score.vocab.ensinado.pos", c("grupo","zona.participante"), aov, pwcs, covar = "score.vocab.ensinado.pre",
    theme = "classic", color = color[["grupo:zona.participante"]],
    subtitle = which(aov$Effect == "grupo:zona.participante"))
}
if (length(unique(pdat[["zona.participante"]])) >= 2) {
  plots[["grupo:zona.participante"]] + ggplot2::ylab("Vocabulary taught") +
  ggplot2::scale_x_discrete(labels=c('pre', 'pos')) +
  if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}
## Warning: No shared levels found between `names(values)` of the manual scale and the data's colour values.

if (length(unique(pdat[["zona.participante"]])) >= 2) {
  plots <- twoWayAncovaBoxPlots(
    wdat.long, "score.vocab.ensinado", c("grupo","zona.participante"), aov, pwc.long,
    pre.post = "time",
    theme = "classic", color = color$prepost)
}
if (length(unique(pdat[["zona.participante"]])) >= 2) 
  plots[["grupo:zona.participante"]] + ggplot2::ylab("Vocabulary taught") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)

Checking linearity assumption

if (length(unique(pdat[["zona.participante"]])) >= 2) {
  ggscatter(wdat, x = "score.vocab.ensinado.pre", y = "score.vocab.ensinado.pos", size = 0.5,
            facet.by = c("grupo","zona.participante"), add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"))
    ) + ggplot2::ylab("Vocabulary taught") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["zona.participante"]])) >= 2) {
  ggscatter(wdat, x = "score.vocab.ensinado.pre", y = "score.vocab.ensinado.pos", size = 0.5,
            color = "grupo", facet.by = "zona.participante", add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = grupo)
    ) +
    ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo:zona.participante"))) +
    ggplot2::scale_color_manual(values = color[["grupo"]]) +
    ggplot2::ylab("Vocabulary taught") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["zona.participante"]])) >= 2) {
  ggscatter(wdat, x = "score.vocab.ensinado.pre", y = "score.vocab.ensinado.pos", size = 0.5,
            color = "zona.participante", facet.by = "grupo", add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = zona.participante)
    ) +
    ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo:zona.participante"))) +
    ggplot2::scale_color_manual(values = color[["zona.participante"]]) +
    ggplot2::ylab("Vocabulary taught") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

Checking normality and homogeneity

if (length(unique(pdat[["zona.participante"]])) >= 2) 
  res <- augment(lm(score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo*zona.participante, data = wdat))
if (length(unique(pdat[["zona.participante"]])) >= 2)
  shapiro_test(res$.resid)
## # A tibble: 1 × 3
##   variable   statistic p.value
##   <chr>          <dbl>   <dbl>
## 1 res$.resid     0.969  0.0856
if (length(unique(pdat[["zona.participante"]])) >= 2) 
  levene_test(res, .resid ~ grupo*zona.participante)
## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     3    65     0.764 0.518

ANCOVA and Pairwise for two factors grupo:zona.escola

Without remove non-normal data

pdat = remove_group_data(dat[!is.na(dat[["grupo"]]) & !is.na(dat[["zona.escola"]]),],
                         "score.vocab.ensinado.pos", c("grupo","zona.escola"))
pdat = pdat[pdat[["zona.escola"]] %in% do.call(
  intersect, lapply(unique(pdat[["grupo"]]), FUN = function(x) {
    unique(pdat[["zona.escola"]][which(pdat[["grupo"]] == x)])
  })),]
pdat[["grupo"]] = factor(pdat[["grupo"]], level[["grupo"]])
pdat[["zona.escola"]] = factor(
  pdat[["zona.escola"]],
  level[["zona.escola"]][level[["zona.escola"]] %in% unique(pdat[["zona.escola"]])])

pdat.long <- rbind(pdat[,c("id","grupo","zona.escola")], pdat[,c("id","grupo","zona.escola")])
pdat.long[["time"]] <- c(rep("pre", nrow(pdat)), rep("pos", nrow(pdat)))
pdat.long[["time"]] <- factor(pdat.long[["time"]], c("pre","pos"))
pdat.long[["score.vocab.ensinado"]] <- c(pdat[["score.vocab.ensinado.pre"]], pdat[["score.vocab.ensinado.pos"]])

if (length(unique(pdat[["zona.escola"]])) >= 2) {
  aov = anova_test(pdat, score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo*zona.escola)
  laov[["grupo:zona.escola"]] <- get_anova_table(aov)
}
if (length(unique(pdat[["zona.escola"]])) >= 2) {
  pwcs <- list()
  pwcs[["zona.escola"]] <- emmeans_test(
    group_by(pdat, grupo), score.vocab.ensinado.pos ~ zona.escola,
    covariate = score.vocab.ensinado.pre, p.adjust.method = "bonferroni")
  pwcs[["grupo"]] <- emmeans_test(
    group_by(pdat, zona.escola), score.vocab.ensinado.pos ~ grupo,
    covariate = score.vocab.ensinado.pre, p.adjust.method = "bonferroni")
  
  pwc <- plyr::rbind.fill(pwcs[["grupo"]], pwcs[["zona.escola"]])
  pwc <- pwc[,c("grupo","zona.escola", colnames(pwc)[!colnames(pwc) %in% c("grupo","zona.escola")])]
}
if (length(unique(pdat[["zona.escola"]])) >= 2) {
  pwc.long <- emmeans_test(dplyr::group_by_at(pdat.long, c("grupo","zona.escola")),
                           score.vocab.ensinado ~ time,
                           p.adjust.method = "bonferroni")
  lpwc[["grupo:zona.escola"]] <- plyr::rbind.fill(pwc, pwc.long)
}
if (length(unique(pdat[["zona.escola"]])) >= 2) {
  ds <- get.descriptives(pdat, "score.vocab.ensinado.pos", c("grupo","zona.escola"), covar = "score.vocab.ensinado.pre")
  ds <- merge(ds[ds$variable != "score.vocab.ensinado.pre",],
              ds[ds$variable == "score.vocab.ensinado.pre", !colnames(ds) %in% c("variable")],
              by = c("grupo","zona.escola"), all.x = T, suffixes = c("", ".score.vocab.ensinado.pre"))
  ds <- merge(get_emmeans(pwcs[["grupo"]]), ds,
              by = c("grupo","zona.escola"), suffixes = c(".emms", ""))
  ds <- ds[,c("grupo","zona.escola","n","mean.score.vocab.ensinado.pre","se.score.vocab.ensinado.pre","mean","se",
              "emmean","se.emms","conf.low","conf.high")]
  
  colnames(ds) <- c("grupo","zona.escola", "N", paste0(c("M","SE")," (pre)"),
                    paste0(c("M","SE"), " (unadj)"),
                    paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")
  
  lemms[["grupo:zona.escola"]] <- ds
}

Computing ANCOVA and PairWise After removing non-normal data (OK)

if (length(unique(pdat[["zona.escola"]])) >= 2) {
  wdat = pdat 
  
  res = residuals(lm(score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo*zona.escola, data = wdat))
  non.normal = getNonNormal(res, wdat$id, plimit = 0.05)
  
  wdat = wdat[!wdat$id %in% non.normal,]
  
  wdat.long <- rbind(wdat[,c("id","grupo","zona.escola")], wdat[,c("id","grupo","zona.escola")])
  wdat.long[["time"]] <- c(rep("pre", nrow(wdat)), rep("pos", nrow(wdat)))
  wdat.long[["time"]] <- factor(wdat.long[["time"]], c("pre","pos"))
  wdat.long[["score.vocab.ensinado"]] <- c(wdat[["score.vocab.ensinado.pre"]], wdat[["score.vocab.ensinado.pos"]])
  
  
  ldat[["grupo:zona.escola"]] = wdat
  
  (non.normal)
}
## [1] "P264" "P286" "P78"  "P290" "P85"
if (length(unique(pdat[["zona.escola"]])) >= 2) {
  aov = anova_test(wdat, score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo*zona.escola)
  laov[["grupo:zona.escola"]] <- merge(get_anova_table(aov), laov[["grupo:zona.escola"]],
                                         by="Effect", suffixes = c("","'"))
  df = get_anova_table(aov)
}
Effect DFn DFd F p p<.05 ges
score.vocab.ensinado.pre 1 82 1.572 0.213 0.019
grupo 1 82 2.241 0.138 0.027
zona.escola 1 82 0.019 0.890 0.000
grupo:zona.escola 1 82 0.307 0.581 0.004
if (length(unique(pdat[["zona.escola"]])) >= 2) {
  pwcs <- list()
  pwcs[["zona.escola"]] <- emmeans_test(
    group_by(wdat, grupo), score.vocab.ensinado.pos ~ zona.escola,
    covariate = score.vocab.ensinado.pre, p.adjust.method = "bonferroni")
  pwcs[["grupo"]] <- emmeans_test(
    group_by(wdat, zona.escola), score.vocab.ensinado.pos ~ grupo,
    covariate = score.vocab.ensinado.pre, p.adjust.method = "bonferroni")
  
  pwc <- plyr::rbind.fill(pwcs[["grupo"]], pwcs[["zona.escola"]])
  pwc <- pwc[,c("grupo","zona.escola", colnames(pwc)[!colnames(pwc) %in% c("grupo","zona.escola")])]
}
grupo zona.escola term .y. group1 group2 df statistic p p.adj p.adj.signif
Rural score.vocab.ensinado.pre*grupo score.vocab.ensinado.pos Controle Experimental 82 -0.343 0.733 0.733 ns
Urbana score.vocab.ensinado.pre*grupo score.vocab.ensinado.pos Controle Experimental 82 -1.561 0.122 0.122 ns
Controle score.vocab.ensinado.pre*zona.escola score.vocab.ensinado.pos Rural Urbana 82 0.290 0.773 0.773 ns
Experimental score.vocab.ensinado.pre*zona.escola score.vocab.ensinado.pos Rural Urbana 82 -0.493 0.623 0.623 ns
if (length(unique(pdat[["zona.escola"]])) >= 2) {
  pwc.long <- emmeans_test(dplyr::group_by_at(wdat.long, c("grupo","zona.escola")),
                           score.vocab.ensinado ~ time,
                           p.adjust.method = "bonferroni")
  lpwc[["grupo:zona.escola"]] <- merge(plyr::rbind.fill(pwc, pwc.long),
                                         lpwc[["grupo:zona.escola"]],
                                         by=c("grupo","zona.escola","term",".y.","group1","group2"),
                                         suffixes = c("","'"))
}
grupo zona.escola term .y. group1 group2 df statistic p p.adj p.adj.signif
Controle Rural time score.vocab.ensinado pre pos 166 0.000 1.000 1.000 ns
Controle Urbana time score.vocab.ensinado pre pos 166 1.014 0.312 0.312 ns
Experimental Rural time score.vocab.ensinado pre pos 166 0.519 0.604 0.604 ns
Experimental Urbana time score.vocab.ensinado pre pos 166 -0.184 0.855 0.855 ns
if (length(unique(pdat[["zona.escola"]])) >= 2) {
  ds <- get.descriptives(wdat, "score.vocab.ensinado.pos", c("grupo","zona.escola"), covar = "score.vocab.ensinado.pre")
  ds <- merge(ds[ds$variable != "score.vocab.ensinado.pre",],
              ds[ds$variable == "score.vocab.ensinado.pre", !colnames(ds) %in% c("variable")],
              by = c("grupo","zona.escola"), all.x = T, suffixes = c("", ".score.vocab.ensinado.pre"))
  ds <- merge(get_emmeans(pwcs[["grupo"]]), ds,
              by = c("grupo","zona.escola"), suffixes = c(".emms", ""))
  ds <- ds[,c("grupo","zona.escola","n","mean.score.vocab.ensinado.pre","se.score.vocab.ensinado.pre",
              "mean","se","emmean","se.emms","conf.low","conf.high")]
  
  colnames(ds) <- c("grupo","zona.escola", "N", paste0(c("M","SE")," (pre)"),
                    paste0(c("M","SE"), " (unadj)"),
                    paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")
  
  lemms[["grupo:zona.escola"]] <- merge(ds, lemms[["grupo:zona.escola"]],
                                          by=c("grupo","zona.escola"), suffixes = c("","'"))
}
grupo zona.escola N M (pre) SE (pre) M (unadj) SE (unadj) M (adj) SE (adj) conf.low conf.high
Controle Rural 12 2.250 0.463 2.250 0.218 2.294 0.435 1.429 3.159
Controle Urbana 36 2.500 0.257 2.139 0.268 2.149 0.250 1.651 2.647
Experimental Rural 13 2.846 0.317 2.538 0.538 2.501 0.417 1.671 3.332
Experimental Urbana 26 2.692 0.313 2.769 0.268 2.753 0.295 2.167 3.340

Plots for ancova

if (length(unique(pdat[["zona.escola"]])) >= 2) {
  ggPlotAoC2(pwcs, "grupo", "zona.escola", aov, ylab = "Vocabulary taught",
             subtitle = which(aov$Effect == "grupo:zona.escola"), addParam = "errorbar") +
    ggplot2::scale_color_manual(values = color[["zona.escola"]]) +
    ggplot2::ylab("Vocabulary taught") +
    if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

if (length(unique(pdat[["zona.escola"]])) >= 2) {
  ggPlotAoC2(pwcs, "zona.escola", "grupo", aov, ylab = "Vocabulary taught",
               subtitle = which(aov$Effect == "grupo:zona.escola"), addParam = "errorbar") +
    ggplot2::scale_color_manual(values = color[["grupo"]]) +
    ggplot2::ylab("Vocabulary taught") +
    if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

if (length(unique(pdat[["zona.escola"]])) >= 2) {
  plots <- twoWayAncovaBoxPlots(
    wdat, "score.vocab.ensinado.pos", c("grupo","zona.escola"), aov, pwcs, covar = "score.vocab.ensinado.pre",
    theme = "classic", color = color[["grupo:zona.escola"]],
    subtitle = which(aov$Effect == "grupo:zona.escola"))
}
if (length(unique(pdat[["zona.escola"]])) >= 2) {
  plots[["grupo:zona.escola"]] + ggplot2::ylab("Vocabulary taught") +
  ggplot2::scale_x_discrete(labels=c('pre', 'pos')) +
  if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}
## Warning: No shared levels found between `names(values)` of the manual scale and the data's colour values.

if (length(unique(pdat[["zona.escola"]])) >= 2) {
  plots <- twoWayAncovaBoxPlots(
    wdat.long, "score.vocab.ensinado", c("grupo","zona.escola"), aov, pwc.long,
    pre.post = "time",
    theme = "classic", color = color$prepost)
}
if (length(unique(pdat[["zona.escola"]])) >= 2) 
  plots[["grupo:zona.escola"]] + ggplot2::ylab("Vocabulary taught") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)

Checking linearity assumption

if (length(unique(pdat[["zona.escola"]])) >= 2) {
  ggscatter(wdat, x = "score.vocab.ensinado.pre", y = "score.vocab.ensinado.pos", size = 0.5,
            facet.by = c("grupo","zona.escola"), add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"))
    ) + ggplot2::ylab("Vocabulary taught") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["zona.escola"]])) >= 2) {
  ggscatter(wdat, x = "score.vocab.ensinado.pre", y = "score.vocab.ensinado.pos", size = 0.5,
            color = "grupo", facet.by = "zona.escola", add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = grupo)
    ) +
    ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo:zona.escola"))) +
    ggplot2::scale_color_manual(values = color[["grupo"]]) +
    ggplot2::ylab("Vocabulary taught") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["zona.escola"]])) >= 2) {
  ggscatter(wdat, x = "score.vocab.ensinado.pre", y = "score.vocab.ensinado.pos", size = 0.5,
            color = "zona.escola", facet.by = "grupo", add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = zona.escola)
    ) +
    ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo:zona.escola"))) +
    ggplot2::scale_color_manual(values = color[["zona.escola"]]) +
    ggplot2::ylab("Vocabulary taught") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

Checking normality and homogeneity

if (length(unique(pdat[["zona.escola"]])) >= 2) 
  res <- augment(lm(score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo*zona.escola, data = wdat))
if (length(unique(pdat[["zona.escola"]])) >= 2)
  shapiro_test(res$.resid)
## # A tibble: 1 × 3
##   variable   statistic p.value
##   <chr>          <dbl>   <dbl>
## 1 res$.resid     0.984   0.365
if (length(unique(pdat[["zona.escola"]])) >= 2) 
  levene_test(res, .resid ~ grupo*zona.escola)
## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     3    83      2.00 0.120

ANCOVA and Pairwise for two factors grupo:score.vocab.ensinado.quintile

Without remove non-normal data

pdat = remove_group_data(dat[!is.na(dat[["grupo"]]) & !is.na(dat[["score.vocab.ensinado.quintile"]]),],
                         "score.vocab.ensinado.pos", c("grupo","score.vocab.ensinado.quintile"))
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning:
## ! There was 1 warning in `mutate()`.
## ℹ In argument: `ci = abs(stats::qt(alpha/2, .data$n - 1) * .data$se)`.
## Caused by warning in `stats::qt()`:
## ! NaNs produced
pdat = pdat[pdat[["score.vocab.ensinado.quintile"]] %in% do.call(
  intersect, lapply(unique(pdat[["grupo"]]), FUN = function(x) {
    unique(pdat[["score.vocab.ensinado.quintile"]][which(pdat[["grupo"]] == x)])
  })),]
pdat[["grupo"]] = factor(pdat[["grupo"]], level[["grupo"]])
pdat[["score.vocab.ensinado.quintile"]] = factor(
  pdat[["score.vocab.ensinado.quintile"]],
  level[["score.vocab.ensinado.quintile"]][level[["score.vocab.ensinado.quintile"]] %in% unique(pdat[["score.vocab.ensinado.quintile"]])])

pdat.long <- rbind(pdat[,c("id","grupo","score.vocab.ensinado.quintile")], pdat[,c("id","grupo","score.vocab.ensinado.quintile")])
pdat.long[["time"]] <- c(rep("pre", nrow(pdat)), rep("pos", nrow(pdat)))
pdat.long[["time"]] <- factor(pdat.long[["time"]], c("pre","pos"))
pdat.long[["score.vocab.ensinado"]] <- c(pdat[["score.vocab.ensinado.pre"]], pdat[["score.vocab.ensinado.pos"]])

if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) {
  aov = anova_test(pdat, score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo*score.vocab.ensinado.quintile)
  laov[["grupo:score.vocab.ensinado.quintile"]] <- get_anova_table(aov)
}
if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) {
  pwcs <- list()
  pwcs[["score.vocab.ensinado.quintile"]] <- emmeans_test(
    group_by(pdat, grupo), score.vocab.ensinado.pos ~ score.vocab.ensinado.quintile,
    covariate = score.vocab.ensinado.pre, p.adjust.method = "bonferroni")
  pwcs[["grupo"]] <- emmeans_test(
    group_by(pdat, score.vocab.ensinado.quintile), score.vocab.ensinado.pos ~ grupo,
    covariate = score.vocab.ensinado.pre, p.adjust.method = "bonferroni")
  
  pwc <- plyr::rbind.fill(pwcs[["grupo"]], pwcs[["score.vocab.ensinado.quintile"]])
  pwc <- pwc[,c("grupo","score.vocab.ensinado.quintile", colnames(pwc)[!colnames(pwc) %in% c("grupo","score.vocab.ensinado.quintile")])]
}
if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) {
  pwc.long <- emmeans_test(dplyr::group_by_at(pdat.long, c("grupo","score.vocab.ensinado.quintile")),
                           score.vocab.ensinado ~ time,
                           p.adjust.method = "bonferroni")
  lpwc[["grupo:score.vocab.ensinado.quintile"]] <- plyr::rbind.fill(pwc, pwc.long)
}
if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) {
  ds <- get.descriptives(pdat, "score.vocab.ensinado.pos", c("grupo","score.vocab.ensinado.quintile"), covar = "score.vocab.ensinado.pre")
  ds <- merge(ds[ds$variable != "score.vocab.ensinado.pre",],
              ds[ds$variable == "score.vocab.ensinado.pre", !colnames(ds) %in% c("variable")],
              by = c("grupo","score.vocab.ensinado.quintile"), all.x = T, suffixes = c("", ".score.vocab.ensinado.pre"))
  ds <- merge(get_emmeans(pwcs[["grupo"]]), ds,
              by = c("grupo","score.vocab.ensinado.quintile"), suffixes = c(".emms", ""))
  ds <- ds[,c("grupo","score.vocab.ensinado.quintile","n","mean.score.vocab.ensinado.pre","se.score.vocab.ensinado.pre","mean","se",
              "emmean","se.emms","conf.low","conf.high")]
  
  colnames(ds) <- c("grupo","score.vocab.ensinado.quintile", "N", paste0(c("M","SE")," (pre)"),
                    paste0(c("M","SE"), " (unadj)"),
                    paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")
  
  lemms[["grupo:score.vocab.ensinado.quintile"]] <- ds
}

Computing ANCOVA and PairWise After removing non-normal data (OK)

if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) {
  wdat = pdat 
  
  res = residuals(lm(score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo*score.vocab.ensinado.quintile, data = wdat))
  non.normal = getNonNormal(res, wdat$id, plimit = 0.05)
  
  wdat = wdat[!wdat$id %in% non.normal,]
  
  wdat.long <- rbind(wdat[,c("id","grupo","score.vocab.ensinado.quintile")], wdat[,c("id","grupo","score.vocab.ensinado.quintile")])
  wdat.long[["time"]] <- c(rep("pre", nrow(wdat)), rep("pos", nrow(wdat)))
  wdat.long[["time"]] <- factor(wdat.long[["time"]], c("pre","pos"))
  wdat.long[["score.vocab.ensinado"]] <- c(wdat[["score.vocab.ensinado.pre"]], wdat[["score.vocab.ensinado.pos"]])
  
  
  ldat[["grupo:score.vocab.ensinado.quintile"]] = wdat
  
  (non.normal)
}
##  [1] "P255" "P286" "P300" "P264" "P76"  "P69"  "P66"  "P294" "P274" "P262" "P167" "P210" "P73"  "P296" "P78"  "P81" 
## [17] "P97"  "P272" "P80"  "P281" "P83"  "P270" "P261" "P256" "P271" "P283" "P236" "P64"  "P277" "P245" "P74"  "P258"
## [33] "P234" "P250" "P93"  "P290" "P251" "P230" "P284" "P273" "P303" "P171" "P297" "P211" "P85"  "P92"  "P65"  "P288"
if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) {
  aov = anova_test(wdat, score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo*score.vocab.ensinado.quintile)
  laov[["grupo:score.vocab.ensinado.quintile"]] <- merge(get_anova_table(aov), laov[["grupo:score.vocab.ensinado.quintile"]],
                                         by="Effect", suffixes = c("","'"))
  df = get_anova_table(aov)
}
Effect DFn DFd F p p<.05 ges
score.vocab.ensinado.pre 1 36 3.033 0.090 0.078
grupo 1 36 3.958 0.054 0.099
score.vocab.ensinado.quintile 2 36 2.496 0.097 0.122
grupo:score.vocab.ensinado.quintile 2 36 9.983 0.000 * 0.357
if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) {
  pwcs <- list()
  pwcs[["score.vocab.ensinado.quintile"]] <- emmeans_test(
    group_by(wdat, grupo), score.vocab.ensinado.pos ~ score.vocab.ensinado.quintile,
    covariate = score.vocab.ensinado.pre, p.adjust.method = "bonferroni")
  pwcs[["grupo"]] <- emmeans_test(
    group_by(wdat, score.vocab.ensinado.quintile), score.vocab.ensinado.pos ~ grupo,
    covariate = score.vocab.ensinado.pre, p.adjust.method = "bonferroni")
  
  pwc <- plyr::rbind.fill(pwcs[["grupo"]], pwcs[["score.vocab.ensinado.quintile"]])
  pwc <- pwc[,c("grupo","score.vocab.ensinado.quintile", colnames(pwc)[!colnames(pwc) %in% c("grupo","score.vocab.ensinado.quintile")])]
}
grupo score.vocab.ensinado.quintile term .y. group1 group2 df statistic p p.adj p.adj.signif
1st quintile score.vocab.ensinado.pre*grupo score.vocab.ensinado.pos Controle Experimental 36 -3.531 0.001 0.001 **
2nd quintile score.vocab.ensinado.pre*grupo score.vocab.ensinado.pos Controle Experimental 36 -2.353 0.024 0.024 *
3rd quintile score.vocab.ensinado.pre*grupo score.vocab.ensinado.pos Controle Experimental 36 2.455 0.019 0.019 *
Controle score.vocab.ensinado.pre*score.vocab.ensinado.quintile score.vocab.ensinado.pos 1st quintile 2nd quintile 36 0.503 0.618 1.000 ns
Controle score.vocab.ensinado.pre*score.vocab.ensinado.quintile score.vocab.ensinado.pos 1st quintile 3rd quintile 36 0.286 0.777 1.000 ns
Controle score.vocab.ensinado.pre*score.vocab.ensinado.quintile score.vocab.ensinado.pos 2nd quintile 3rd quintile 36 -0.014 0.989 1.000 ns
Experimental score.vocab.ensinado.pre*score.vocab.ensinado.quintile score.vocab.ensinado.pos 1st quintile 2nd quintile 36 1.862 0.071 0.213 ns
Experimental score.vocab.ensinado.pre*score.vocab.ensinado.quintile score.vocab.ensinado.pos 1st quintile 3rd quintile 36 3.301 0.002 0.007 **
Experimental score.vocab.ensinado.pre*score.vocab.ensinado.quintile score.vocab.ensinado.pos 2nd quintile 3rd quintile 36 3.490 0.001 0.004 **
if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) {
  pwc.long <- emmeans_test(dplyr::group_by_at(wdat.long, c("grupo","score.vocab.ensinado.quintile")),
                           score.vocab.ensinado ~ time,
                           p.adjust.method = "bonferroni")
  lpwc[["grupo:score.vocab.ensinado.quintile"]] <- merge(plyr::rbind.fill(pwc, pwc.long),
                                         lpwc[["grupo:score.vocab.ensinado.quintile"]],
                                         by=c("grupo","score.vocab.ensinado.quintile","term",".y.","group1","group2"),
                                         suffixes = c("","'"))
}
grupo score.vocab.ensinado.quintile term .y. group1 group2 df statistic p p.adj p.adj.signif
Controle 1st quintile time score.vocab.ensinado pre pos 74 -5.697 0.000 0.000 ****
Controle 2nd quintile time score.vocab.ensinado pre pos 74 -1.415 0.161 0.161 ns
Controle 3rd quintile time score.vocab.ensinado pre pos 74 1.513 0.135 0.135 ns
Experimental 1st quintile time score.vocab.ensinado pre pos 74 -6.356 0.000 0.000 ****
Experimental 2nd quintile time score.vocab.ensinado pre pos 74 -4.431 0.000 0.000 ****
Experimental 3rd quintile time score.vocab.ensinado pre pos 74 3.581 0.001 0.001 ***
if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) {
  ds <- get.descriptives(wdat, "score.vocab.ensinado.pos", c("grupo","score.vocab.ensinado.quintile"), covar = "score.vocab.ensinado.pre")
  ds <- merge(ds[ds$variable != "score.vocab.ensinado.pre",],
              ds[ds$variable == "score.vocab.ensinado.pre", !colnames(ds) %in% c("variable")],
              by = c("grupo","score.vocab.ensinado.quintile"), all.x = T, suffixes = c("", ".score.vocab.ensinado.pre"))
  ds <- merge(get_emmeans(pwcs[["grupo"]]), ds,
              by = c("grupo","score.vocab.ensinado.quintile"), suffixes = c(".emms", ""))
  ds <- ds[,c("grupo","score.vocab.ensinado.quintile","n","mean.score.vocab.ensinado.pre","se.score.vocab.ensinado.pre",
              "mean","se","emmean","se.emms","conf.low","conf.high")]
  
  colnames(ds) <- c("grupo","score.vocab.ensinado.quintile", "N", paste0(c("M","SE")," (pre)"),
                    paste0(c("M","SE"), " (unadj)"),
                    paste0(c("M", "SE"), " (adj)"), "conf.low", "conf.high")
  
  lemms[["grupo:score.vocab.ensinado.quintile"]] <- merge(ds, lemms[["grupo:score.vocab.ensinado.quintile"]],
                                          by=c("grupo","score.vocab.ensinado.quintile"), suffixes = c("","'"))
}
grupo score.vocab.ensinado.quintile N M (pre) SE (pre) M (unadj) SE (unadj) M (adj) SE (adj) conf.low conf.high
Controle 1st quintile 10 0.600 0.163 2.400 0.221 3.160 0.507 2.131 4.189
Controle 2nd quintile 8 2.375 0.183 2.875 0.227 2.864 0.290 2.276 3.451
Controle 3rd quintile 7 4.286 0.286 3.714 0.565 2.873 0.574 1.709 4.037
Experimental 1st quintile 3 0.667 0.333 4.333 0.333 5.064 0.632 3.782 6.346
Experimental 2nd quintile 10 2.400 0.163 3.800 0.200 3.778 0.259 3.252 4.304
Experimental 3rd quintile 5 4.000 0.000 2.400 0.245 1.683 0.551 0.565 2.800

Plots for ancova

if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) {
  ggPlotAoC2(pwcs, "grupo", "score.vocab.ensinado.quintile", aov, ylab = "Vocabulary taught",
             subtitle = which(aov$Effect == "grupo:score.vocab.ensinado.quintile"), addParam = "errorbar") +
    ggplot2::scale_color_manual(values = color[["score.vocab.ensinado.quintile"]]) +
    ggplot2::ylab("Vocabulary taught") +
    if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) {
  ggPlotAoC2(pwcs, "score.vocab.ensinado.quintile", "grupo", aov, ylab = "Vocabulary taught",
               subtitle = which(aov$Effect == "grupo:score.vocab.ensinado.quintile"), addParam = "errorbar") +
    ggplot2::scale_color_manual(values = color[["grupo"]]) +
    ggplot2::ylab("Vocabulary taught") +
    if (ymin.ci < ymax.ci) ggplot2::ylim(ymin.ci, ymax.ci)
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) {
  plots <- twoWayAncovaBoxPlots(
    wdat, "score.vocab.ensinado.pos", c("grupo","score.vocab.ensinado.quintile"), aov, pwcs, covar = "score.vocab.ensinado.pre",
    theme = "classic", color = color[["grupo:score.vocab.ensinado.quintile"]],
    subtitle = which(aov$Effect == "grupo:score.vocab.ensinado.quintile"))
}
if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) {
  plots[["grupo:score.vocab.ensinado.quintile"]] + ggplot2::ylab("Vocabulary taught") +
  ggplot2::scale_x_discrete(labels=c('pre', 'pos')) +
  if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}
## Warning: No shared levels found between `names(values)` of the manual scale and the data's colour values.

if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) {
  plots <- twoWayAncovaBoxPlots(
    wdat.long, "score.vocab.ensinado", c("grupo","score.vocab.ensinado.quintile"), aov, pwc.long,
    pre.post = "time",
    theme = "classic", color = color$prepost)
}
if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) 
  plots[["grupo:score.vocab.ensinado.quintile"]] + ggplot2::ylab("Vocabulary taught") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)

Checking linearity assumption

if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) {
  ggscatter(wdat, x = "score.vocab.ensinado.pre", y = "score.vocab.ensinado.pos", size = 0.5,
            facet.by = c("grupo","score.vocab.ensinado.quintile"), add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"))
    ) + ggplot2::ylab("Vocabulary taught") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) {
  ggscatter(wdat, x = "score.vocab.ensinado.pre", y = "score.vocab.ensinado.pos", size = 0.5,
            color = "grupo", facet.by = "score.vocab.ensinado.quintile", add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = grupo)
    ) +
    ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo:score.vocab.ensinado.quintile"))) +
    ggplot2::scale_color_manual(values = color[["grupo"]]) +
    ggplot2::ylab("Vocabulary taught") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) {
  ggscatter(wdat, x = "score.vocab.ensinado.pre", y = "score.vocab.ensinado.pos", size = 0.5,
            color = "score.vocab.ensinado.quintile", facet.by = "grupo", add = "reg.line")+
    stat_regline_equation(
      aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = score.vocab.ensinado.quintile)
    ) +
    ggplot2::labs(subtitle = rstatix::get_test_label(aov, detailed = T, row = which(aov$Effect == "grupo:score.vocab.ensinado.quintile"))) +
    ggplot2::scale_color_manual(values = color[["score.vocab.ensinado.quintile"]]) +
    ggplot2::ylab("Vocabulary taught") +
    if (ymin < ymax) ggplot2::ylim(ymin, ymax)
}

Checking normality and homogeneity

if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) 
  res <- augment(lm(score.vocab.ensinado.pos ~ score.vocab.ensinado.pre + grupo*score.vocab.ensinado.quintile, data = wdat))
if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2)
  shapiro_test(res$.resid)
## # A tibble: 1 × 3
##   variable   statistic p.value
##   <chr>          <dbl>   <dbl>
## 1 res$.resid     0.964   0.202
if (length(unique(pdat[["score.vocab.ensinado.quintile"]])) >= 2) 
  levene_test(res, .resid ~ grupo*score.vocab.ensinado.quintile)
## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     5    37      1.29 0.290

Summary of Results

Descriptive Statistics

df <- get.descriptives(ldat[["grupo"]], c(dv.pre, dv.pos), c("grupo"), 
                       include.global = T, symmetry.test = T, normality.test = F)
df <- plyr::rbind.fill(
  df, do.call(plyr::rbind.fill, lapply(lfatores2, FUN = function(f) {
    if (nrow(dat) > 0 && sum(!is.na(unique(dat[[f]]))) > 1 && paste0("grupo:",f) %in% names(ldat))
      get.descriptives(ldat[[paste0("grupo:",f)]], c(dv.pre,dv.pos), c("grupo", f),
                       symmetry.test = T, normality.test = F)
    }))
)
df <- df[,c(fatores1[fatores1 %in% colnames(df)],"variable",
             colnames(df)[!colnames(df) %in% c(fatores1,"variable")])]
grupo genero zona.participante zona.escola score.vocab.ensinado.quintile variable n mean median min max sd se ci iqr symmetry skewness kurtosis
Controle score.vocab.ensinado.pre 48 2.333 2.0 0 6 1.464 0.211 0.425 2.25 YES 0.183 -0.759
Experimental score.vocab.ensinado.pre 36 2.667 2.0 0 7 1.474 0.246 0.499 2.00 NO 0.619 0.205
score.vocab.ensinado.pre 84 2.476 2.0 0 7 1.468 0.160 0.319 3.00 YES 0.373 -0.170
Controle score.vocab.ensinado.pos 48 2.167 2.0 0 6 1.434 0.207 0.416 2.00 NO 0.646 0.412
Experimental score.vocab.ensinado.pos 36 2.750 2.5 0 5 1.381 0.230 0.467 2.00 YES 0.059 -1.156
score.vocab.ensinado.pos 84 2.417 2.0 0 6 1.433 0.156 0.311 1.25 YES 0.372 -0.369
Controle F score.vocab.ensinado.pre 24 2.625 3.0 0 6 1.527 0.312 0.645 2.25 YES 0.124 -0.676
Controle M score.vocab.ensinado.pre 24 2.042 2.0 0 4 1.367 0.279 0.577 2.00 YES 0.125 -1.360
Experimental F score.vocab.ensinado.pre 16 2.312 2.0 0 5 1.352 0.338 0.721 2.00 YES 0.222 -0.949
Experimental M score.vocab.ensinado.pre 20 2.950 2.5 1 7 1.538 0.344 0.720 2.00 NO 0.738 0.065
Controle F score.vocab.ensinado.pos 24 2.167 2.0 0 6 1.494 0.305 0.631 1.25 NO 1.152 1.165
Controle M score.vocab.ensinado.pos 24 2.167 2.0 0 5 1.404 0.287 0.593 2.00 YES -0.012 -0.909
Experimental F score.vocab.ensinado.pos 16 2.688 2.5 0 5 1.537 0.384 0.819 2.25 YES -0.022 -1.360
Experimental M score.vocab.ensinado.pos 20 2.800 2.5 1 5 1.281 0.287 0.600 2.00 YES 0.211 -1.302
Controle Rural score.vocab.ensinado.pre 12 2.000 2.5 0 4 1.537 0.444 0.977 2.25 YES -0.138 -1.717
Controle Urbana score.vocab.ensinado.pre 26 2.500 2.0 0 6 1.631 0.320 0.659 3.00 YES 0.479 -0.610
Experimental Rural score.vocab.ensinado.pre 15 2.800 3.0 1 5 1.265 0.327 0.700 1.50 YES 0.348 -1.072
Experimental Urbana score.vocab.ensinado.pre 16 2.938 3.0 0 7 1.611 0.403 0.859 2.00 NO 0.543 0.443
Controle Rural score.vocab.ensinado.pos 12 2.833 2.0 1 8 1.801 0.520 1.144 1.00 NO 1.854 2.776
Controle Urbana score.vocab.ensinado.pos 26 2.192 2.0 0 6 1.877 0.368 0.758 2.00 NO 0.780 -0.425
Experimental Rural score.vocab.ensinado.pos 15 2.333 2.0 0 5 1.496 0.386 0.828 2.00 YES 0.062 -1.140
Experimental Urbana score.vocab.ensinado.pos 16 3.000 3.5 0 6 1.826 0.456 0.973 2.25 YES -0.246 -1.234
Controle Rural score.vocab.ensinado.pre 12 2.250 2.0 0 5 1.603 0.463 1.018 2.25 YES 0.235 -1.550
Controle Urbana score.vocab.ensinado.pre 36 2.500 2.0 0 6 1.540 0.257 0.521 2.25 YES 0.251 -0.415
Experimental Rural score.vocab.ensinado.pre 13 2.846 3.0 1 5 1.144 0.317 0.691 2.00 YES 0.272 -1.155
Experimental Urbana score.vocab.ensinado.pre 26 2.692 2.5 0 7 1.594 0.313 0.644 2.75 NO 0.545 -0.045
Controle Rural score.vocab.ensinado.pos 12 2.250 2.0 1 3 0.754 0.218 0.479 1.00 few data 0.000 0.000
Controle Urbana score.vocab.ensinado.pos 36 2.139 2.0 0 6 1.606 0.268 0.543 2.00 NO 0.665 -0.090
Experimental Rural score.vocab.ensinado.pos 13 2.538 2.0 0 6 1.941 0.538 1.173 3.00 YES 0.356 -1.288
Experimental Urbana score.vocab.ensinado.pos 26 2.769 3.0 0 5 1.366 0.268 0.552 2.00 YES -0.142 -1.138
Controle 1st quintile score.vocab.ensinado.pre 10 0.600 1.0 0 1 0.516 0.163 0.369 1.00 few data 0.000 0.000
Controle 2nd quintile score.vocab.ensinado.pre 8 2.375 2.0 2 3 0.518 0.183 0.433 1.00 few data 0.000 0.000
Controle 3rd quintile score.vocab.ensinado.pre 7 4.286 4.0 4 6 0.756 0.286 0.699 0.00 few data 0.000 0.000
Experimental 1st quintile score.vocab.ensinado.pre 3 0.667 1.0 0 1 0.577 0.333 1.434 0.50 few data 0.000 0.000
Experimental 2nd quintile score.vocab.ensinado.pre 10 2.400 2.0 2 3 0.516 0.163 0.369 1.00 few data 0.000 0.000
Experimental 3rd quintile score.vocab.ensinado.pre 5 4.000 4.0 4 4 0.000 0.000 0.000 0.00 few data 0.000 0.000
Controle 1st quintile score.vocab.ensinado.pos 10 2.400 2.0 2 4 0.699 0.221 0.500 0.75 few data 0.000 0.000
Controle 2nd quintile score.vocab.ensinado.pos 8 2.875 3.0 2 4 0.641 0.227 0.536 0.25 few data 0.000 0.000
Controle 3rd quintile score.vocab.ensinado.pos 7 3.714 4.0 2 6 1.496 0.565 1.384 2.00 YES 0.157 -1.643
Experimental 1st quintile score.vocab.ensinado.pos 3 4.333 4.0 4 5 0.577 0.333 1.434 0.50 few data 0.000 0.000
Experimental 2nd quintile score.vocab.ensinado.pos 10 3.800 4.0 3 5 0.632 0.200 0.452 0.75 few data 0.000 0.000
Experimental 3rd quintile score.vocab.ensinado.pos 5 2.400 2.0 2 3 0.548 0.245 0.680 1.00 few data 0.000 0.000

ANCOVA Table Comparison

df <- do.call(plyr::rbind.fill, laov)
df <- df[!duplicated(df$Effect),]
Effect DFn DFd F p p<.05 ges DFn’ DFd’ F’ p’ p<.05’ ges’
1 grupo 1 81 3.174 0.079 0.038 1 89 0.487 0.487 0.005
2 score.vocab.ensinado.pre 1 81 0.429 0.514 0.005 1 89 0.382 0.538 0.004
3 genero 1 79 0.028 0.868 0.000 1 87 0.001 0.975 0.000
5 grupo:genero 1 79 0.002 0.966 0.000 1 87 0.032 0.858 0.000
8 grupo:zona.participante 1 64 2.431 0.124 0.037 1 64 2.431 0.124 0.037
10 zona.participante 1 64 0.022 0.881 0.000 1 64 0.022 0.881 0.000
12 grupo:zona.escola 1 82 0.307 0.581 0.004 1 87 1.523 0.220 0.017
14 zona.escola 1 82 0.019 0.890 0.000 1 87 0.837 0.363 0.010
16 grupo:score.vocab.ensinado.quintile 2 36 9.983 0.000 * 0.357 2 84 2.030 0.138 0.046
18 score.vocab.ensinado.quintile 2 36 2.496 0.097 0.122 2 84 2.586 0.081 0.058

PairWise Table Comparison

df <- do.call(plyr::rbind.fill, lpwc)
df <- df[,c(names(lfatores)[names(lfatores) %in% colnames(df)],
            names(df)[!names(df) %in% c(names(lfatores),"term",".y.")])]
grupo genero zona.participante zona.escola score.vocab.ensinado.quintile group1 group2 df statistic p p.adj p.adj.signif df’ statistic’ p’ p.adj’ p.adj.signif’
Controle pre pos 164 0.567 0.571 0.571 ns 180 0.187 0.852 0.852 ns
Experimental pre pos 164 -0.246 0.806 0.806 ns 180 0.545 0.586 0.586 ns
Controle Experimental 81 -1.782 0.079 0.079 ns 89 -0.698 0.487 0.487 ns
Controle F pre pos 160 1.102 0.272 0.272 ns 176 0.527 0.599 0.599 ns
Controle M pre pos 160 -0.300 0.764 0.764 ns 176 -0.263 0.792 0.792 ns
Controle F M 79 -0.096 0.923 0.923 ns 87 0.099 0.922 0.922 ns
Experimental F pre pos 160 -0.736 0.463 0.463 ns 176 0.000 1.000 1.000 ns
Experimental M pre pos 160 0.329 0.742 0.742 ns 176 0.733 0.465 0.465 ns
Experimental F M 79 -0.141 0.889 0.889 ns 87 -0.155 0.877 0.877 ns
F Controle Experimental 79 -1.169 0.246 0.246 ns 87 -0.348 0.728 0.728 ns
M Controle Experimental 79 -1.280 0.204 0.204 ns 87 -0.618 0.538 0.538 ns
Controle Rural Urbana 64 1.195 0.237 0.237 ns 64 1.195 0.237 0.237 ns
Controle Rural pre pos 130 -1.229 0.221 0.221 ns 130 -1.229 0.221 0.221 ns
Controle Urbana pre pos 130 0.668 0.506 0.506 ns 130 0.668 0.506 0.506 ns
Experimental Rural Urbana 64 -1.008 0.317 0.317 ns 64 -1.008 0.317 0.317 ns
Experimental Rural pre pos 130 0.769 0.443 0.443 ns 130 0.769 0.443 0.443 ns
Experimental Urbana pre pos 130 -0.106 0.915 0.915 ns 130 -0.106 0.915 0.915 ns
Rural Controle Experimental 64 0.951 0.345 0.345 ns 64 0.951 0.345 0.345 ns
Urbana Controle Experimental 64 -1.279 0.206 0.206 ns 64 -1.279 0.206 0.206 ns
Controle Rural Urbana 82 0.290 0.773 0.773 ns 87 1.513 0.134 0.134 ns
Controle Rural pre pos 166 0.000 1.000 1.000 ns 176 -1.177 0.241 0.241 ns
Controle Urbana pre pos 166 1.014 0.312 0.312 ns 176 0.954 0.341 0.341 ns
Experimental Rural Urbana 82 -0.493 0.623 0.623 ns 87 -0.268 0.789 0.789 ns
Experimental Rural pre pos 166 0.519 0.604 0.604 ns 176 0.489 0.626 0.626 ns
Experimental Urbana pre pos 166 -0.184 0.855 0.855 ns 176 0.327 0.744 0.744 ns
Rural Controle Experimental 82 -0.343 0.733 0.733 ns 87 0.672 0.503 0.503 ns
Urbana Controle Experimental 82 -1.561 0.122 0.122 ns 87 -1.222 0.225 0.225 ns
Controle 1st quintile pre pos 74 -5.697 0.000 0.000 **** 170 -2.979 0.003 0.003 **
Controle 2nd quintile pre pos 74 -1.415 0.161 0.161 ns 170 0.248 0.805 0.805 ns
Controle 3rd quintile pre pos 74 1.513 0.135 0.135 ns 170 3.463 0.001 0.001 ***
Controle 1st quintile 2nd quintile 36 0.503 0.618 1.000 ns 84 0.691 0.492 1.000 ns
Controle 1st quintile 3rd quintile 36 0.286 0.777 1.000 ns 84 0.946 0.347 1.000 ns
Controle 2nd quintile 3rd quintile 36 -0.014 0.989 1.000 ns 84 0.872 0.385 1.000 ns
Experimental 1st quintile pre pos 74 -6.356 0.000 0.000 **** 170 -3.010 0.003 0.003 **
Experimental 2nd quintile pre pos 74 -4.431 0.000 0.000 **** 170 -1.823 0.070 0.070 ns
Experimental 3rd quintile pre pos 74 3.581 0.001 0.001 *** 170 5.157 0.000 0.000 ****
Experimental 1st quintile 2nd quintile 36 1.862 0.071 0.213 ns 84 0.500 0.618 1.000 ns
Experimental 1st quintile 3rd quintile 36 3.301 0.002 0.007 ** 84 2.052 0.043 0.130 ns
Experimental 2nd quintile 3rd quintile 36 3.490 0.001 0.004 ** 84 2.772 0.007 0.021 *
1st quintile Controle Experimental 36 -3.531 0.001 0.001 ** 84 -0.893 0.374 0.374 ns
2nd quintile Controle Experimental 36 -2.353 0.024 0.024 * 84 -1.486 0.141 0.141 ns
3rd quintile Controle Experimental 36 2.455 0.019 0.019 * 84 1.262 0.210 0.210 ns

EMMS Table Comparison

df <- do.call(plyr::rbind.fill, lemms)
df[["N-N'"]] <- df[["N"]] - df[["N'"]]
df <- df[,c(names(lfatores)[names(lfatores) %in% colnames(df)],
            names(df)[!names(df) %in% names(lfatores)])]
grupo genero zona.participante zona.escola score.vocab.ensinado.quintile N M (pre) SE (pre) M (unadj) SE (unadj) M (adj) SE (adj) conf.low conf.high N’ M (pre)’ SE (pre)’ M (unadj)’ SE (unadj)’ M (adj)’ SE (adj)’ conf.low’ conf.high’ N-N’
Controle 48 2.333 0.211 2.167 0.207 2.177 0.205 1.769 2.585 50 2.420 0.216 2.360 0.242 2.374 0.243 1.892 2.857 -2
Experimental 36 2.667 0.246 2.750 0.230 2.737 0.237 2.265 3.208 42 2.833 0.223 2.643 0.261 2.626 0.265 2.099 3.153 -6
Controle F 24 2.625 0.312 2.167 0.305 2.156 0.293 1.573 2.740 25 2.640 0.299 2.400 0.374 2.398 0.346 1.711 3.085 -1
Controle M 24 2.042 0.279 2.167 0.287 2.197 0.297 1.606 2.787 25 2.200 0.311 2.320 0.315 2.349 0.349 1.655 3.044 -1
Experimental F 16 2.312 0.338 2.688 0.384 2.699 0.359 1.984 3.413 19 2.579 0.327 2.579 0.421 2.581 0.397 1.793 3.369 -3
Experimental M 20 2.950 0.344 2.800 0.287 2.767 0.325 2.120 3.414 23 3.043 0.305 2.696 0.335 2.665 0.364 1.940 3.389 -3
Controle Rural 12 2.000 0.444 2.833 0.520 2.948 0.516 1.918 3.978 12 2.000 0.444 2.833 0.520 2.948 0.516 1.918 3.978 0
Controle Urbana 26 2.500 0.320 2.192 0.368 2.208 0.346 1.517 2.899 26 2.500 0.320 2.192 0.368 2.208 0.346 1.517 2.899 0
Experimental Rural 15 2.800 0.327 2.333 0.386 2.290 0.456 1.378 3.201 15 2.800 0.327 2.333 0.386 2.290 0.456 1.378 3.201 0
Experimental Urbana 16 2.938 0.403 3.000 0.456 2.929 0.444 2.043 3.816 16 2.938 0.403 3.000 0.456 2.929 0.444 2.043 3.816 0
Controle Rural 12 2.250 0.463 2.250 0.218 2.294 0.435 1.429 3.159 14 2.214 0.408 2.929 0.508 2.963 0.458 2.052 3.874 -2
Controle Urbana 36 2.500 0.257 2.139 0.268 2.149 0.250 1.651 2.647 36 2.500 0.257 2.139 0.268 2.148 0.285 1.583 2.714 0
Experimental Rural 13 2.846 0.317 2.538 0.538 2.501 0.417 1.671 3.332 13 2.846 0.317 2.538 0.538 2.518 0.474 1.576 3.460 0
Experimental Urbana 26 2.692 0.313 2.769 0.268 2.753 0.295 2.167 3.340 29 2.828 0.294 2.690 0.298 2.671 0.318 2.039 3.303 -3
Controle 1st quintile 10 0.600 0.163 2.400 0.221 3.160 0.507 2.131 4.189 16 0.688 0.120 2.000 0.387 3.027 0.765 1.505 4.549 -6
Controle 2nd quintile 8 2.375 0.183 2.875 0.227 2.864 0.290 2.276 3.451 21 2.524 0.112 2.429 0.406 2.449 0.364 1.725 3.172 -13
Controle 3rd quintile 7 4.286 0.286 3.714 0.565 2.873 0.574 1.709 4.037 13 4.385 0.213 2.692 0.472 1.692 0.778 0.145 3.238 -6
Experimental 1st quintile 3 0.667 0.333 4.333 0.333 5.064 0.632 3.782 6.346 8 0.875 0.125 2.750 0.491 3.674 0.825 2.034 5.315 -5
Experimental 2nd quintile 10 2.400 0.163 3.800 0.200 3.778 0.259 3.252 4.304 19 2.421 0.116 3.158 0.377 3.234 0.385 2.468 4.000 -9
Experimental 3rd quintile 5 4.000 0.000 2.400 0.245 1.683 0.551 0.565 2.800 14 4.214 0.114 1.786 0.459 0.879 0.721 -0.555 2.312 -9